www.gusucode.com > 超声波测量以及形成图像 对相关信号进行模拟仿真 > 超声波测量以及形成图像 对相关信号进行模拟仿真/digital holograpy/prog/FAFS.m
function [ varargout ] = FAFS( U,varargin ) %FAFS Fourier transform with Arbitrary Frequency Sampling % Syntax: % [F,dfxy] = FAFS(U,dxy,b,FSN,f0) % [F,...] = FAFS(U,...) % % U : the wavefront on the object plane % size of U must be even % F : frequency distribution of U % dxy : sampling interval of U, dxy=[dx,dy] % dfxy : sampling interval in frequency domain, dfxy=[dfx,dfy] % b : the parameter to set the dfx/dfy, b=[b1,b2] % FSN : the number of sampling points on the frequency plane % FSN=[FSNy,FSNx] <---ATTENTION % f0 : the center of the sampling region on frequency plane % f0=[f0x,f0y] % % if there is no output, frequency image will be displayed % else, no image will be displayed % % the origin of coordinates is at M+1/2,N/2+1 % % the relation between sampling distance on the two planes is: % dfx=1/dx/N/b(1),dfy=1/dy/M/b(2) % error(nargchk(1,5,nargin)) if nargout>2 error('Too many output arguments') end [M,N]=size(U); switch nargin case 1 dxy=[1,1]; b=[1,1]; FSN=[M,N]; f0=[0,0]; case 2 dxy=varargin{1}; b=[1,1]; FSN=[M,N]; f0=[0,0]; case 3 dxy=varargin{1}; b=varargin{2}; FSN=[M,N]; f0=[0,0]; case 4 dxy=varargin{1}; b=varargin{2}; FSN=varargin{3}; f0=[0,0]; case 5 dxy=varargin{1}; b=varargin{2}; FSN=varargin{3}; f0=varargin{4}; end %---------------------------------------- [P,Q]=meshgrid(0:N-1,0:M-1); phase=exp(i*pi*((P+1)*FSN(2)/b(1)/N+(Q+1)*FSN(1)/b(2)/M-2*f0(1)*dxy(1)*(P+1)+2*f0(2)*dxy(2)*(Q+1)-P.^2/b(1)/N-Q.^2/b(2)/M)); U=U.*phase; clear phase; U=paste(zeros(FSN(1)+M-1,FSN(2)+N-1),U,M/2+1,N/2+1); FU=fft2(U); clear U %---------------------------------------- [P,Q]=meshgrid(1-N:FSN(2)-1,1-M:FSN(1)-1); h=exp(i*pi*(P.^2/b(1)/N+Q.^2/b(2)/M)); Fh=fft2(h); clear h %---------------------------------------- FU=FU.*Fh; clear Fh %---------------------------------------- [P,Q]=meshgrid(0:FSN(2)+N-2,0:FSN(1)+M-2); fp=exp(i*2*pi*((N-1)/(FSN(2)+N-1)*P+(M-1)/(FSN(1)+M-1)*Q)); FU=FU.*fp; clear fp P Q %---------------------------------------- F=ifft2(FU); clear FU F=F(1:FSN(1),1:FSN(2)); %---------------------------------------- [P,Q]=meshgrid(1:FSN(2),1:FSN(1)); phase=exp(i*pi*((P-FSN(2)/2-(FSN(2)+N)/N)/b(1)+(Q-FSN(1)/2-(FSN(1)+M)/M)/b(2)+... 2*f0(1)*(N/2+1)*dxy(1)-2*f0(2)*(M/2+1)*dxy(2)-... (P-1).^2/b(1)/N-(Q-1).^2/b(2)/M))*dxy(1)*dxy(2); clear P Q F=F.*phase; %---------------------------------------- switch nargout case 0 f=abs(F); fx=(-FSN(2)/2:FSN(2)/2-1)/dxy(1)/N/b(1)+f0(1); fy=(FSN(1)/2:-1:1-FSN(1)/2)/dxy(2)/M/b(2)+f0(2); imagesc(fx,fy,f);colormap(gray);axis xy image; xlabel('fx');ylabel('fy');title('Amplitude of Frequency Spectrum'); case 1 varargout{1}=F; case 2 varargout{1}=F; varargout{2}=[1/dxy(1)/N/b(1),1/dxy(2)/M/b(2)]; end